Single and two-particle energy gaps across the disorder-driven 
superconductor-insulator transition 



Karim Bouadim, Yen Lee Loh, Mohit Randeria, and Nandini Trivedi 

Department of Physics, The Ohio State University, Columbus, OH 43210, USA 

The competition between superconductivity and localization raises profound questions in con- 
densed matter physics. In spite of decades of research, the mechanism of the superconductor- 
insulator transition (SIT) and the nature of the insulator are not understood. We use quantum 
Monte Carlo simulations that treat, on an equal footing, inhomogeneous amplitude variations and 
phase fluctuations, a major advance over previous theories. We gain new microscopic insights and 
make testable predictions for local spectroscopic probes. The energy gap in the density of states sur- 
vives across the transition, but coherence peaks exist only in the superconductor. A characteristic 
pseudogap persists above the critical disorder and critical temperature, in contrast to conventional 
theories. Surprisingly, the insulator has a two-particle gap scale that vanishes at the SIT, despite a 
robust single-particle gap. 
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Attractive interactions between electrons lead to su- 
perconductivity, a spectacular example of long range or- 
der in physics, while disorder leads to localization of elec- 
tronic states. One of the most fascinating examples of the 
interplay between the effects of interactions and disorder, 
is the destruction of superconductivity in thin films with 
increasing disorder (or decreasing thickness or increas- 
ing magnetic field), and the resulting superconductor-to- 
insulator transition (SIT) [1-9]. 

It was recognized long ago that s-wave superconductiv- 
ity is remarkably robust against weak disorder [10, 11]. 
But strong disorder, where the SIT occurs, is very hard 
to treat theoretically in an interacting system. Early 
work on critical phenomena [12] at the SIT used bosonic 
models to describe fluctuations of the phase of the su- 
perconducting order parameter. A more microscopic ap- 
proach, involving the fermionic degrees of freedom, leads 
to a highly inhomogeneous state [13-16]. Even though 
the disorder potential varies on an atomic scale, the su- 
perconducting pairing amplitude forms "self-organized" 
puddles on the scale of the coherence length and results in 
large suppression of the superfluid density [13, 14]. How- 
ever, the inhomogeneous mean field theory by itself fails 
to capture the SIT, which necessarily requires inclusion 
of phase fluctuations. 

In this paper we make a major advance by using quan- 
tum Monte Carlo (QMC) simulations of a microscopic 
model for an s-wave superconductor (SC) in a random 
potential, which treats on an equal footing the inhomo- 
geneous variations of the pairing amplitude and thermal 
and quantum phase fluctuations. What is new in our pa- 
per are calculations of one-particle and two-particle spec- 
tral functions across the SIT which allows us to answer 
the questions about the mechanism of the transition, the 
evolution of the local excitation spectrum across the SIT 
and the nature of the resulting insulator. 

Our main results are as follows: 
(1) Single-particle gap: At T = the gap in the single- 



particle density of states (DOS) survives through the 
SIT, so that one goes from a gapped superconductor to a 
gapped insulator. Although both states are highly inho- 
mogeneous, the local density of states (LDOS) is gapped 
at every site. 

(2) Coherence peaks: Coherence peaks - characteristic 
pile-ups in the DOS at the gap edges - are strongly corre- 
lated with superconducting order and vanish as the tem- 
perature is raised above T c , or as the disorder is increased 
across the SIT. 

(3) Pseudogap: Near the SIT, a pseudogap - a suppres- 
sion in the low-energy DOS - persists well above the su- 
perconducting T c , in marked deviation from BCS theory. 
A pseudogap also exists at finite temperatures in the in- 
sulating state. 

(4) Two-particle gap: There is a characteristic energy 
scale c^pair to insert a pair in the insulator that collapses 
upon approaching the SIT from the insulating side. How- 
ever, the two-particle spectral function does not have a 
hard gap since rare regions lead to a very small spectral 
weight at low energies. 

Our results thus provide a complete description of the 
phases and the quantum critical region bounded by T c 
in the superconductor and cj pa i r in the insulator. Fur- 
thermore, our predictions for the local tunneling density 
of states and the dynamical pair susceptibility as a func- 
tion of temperature and disorder have the potential to 
guide future experiments using scanning tunneling spec- 
troscopy (STS) [17-20] and other dynamical probes [21]. 

Model and methods: To model the competition be- 
tween superconductivity and localization that leads to 
the SIT in quench-condensed films with thicknesses less 
than the coherence length, we take the simplest lattice 
Hamiltonian that has an s-wave superconducting ground 
state in the absence of disorder (V = 0) and exhibits 
Anderson localization when the attractive interaction is 
turned off (U = 0). This is the two-dimensional attrac- 
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tive Hubbard model in a random potential: 
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with lattice sites R and R/, spin indices a =t or ^, 
fermion creation and annihilation operators and c , 
number operators nR a = Cp^Cp^, hopping t between 
neighboring sites (RR'}, and chemical potential fi. Vr is 
a random potential at each site drawn from the uniform 
distribution on [— V, +V], and \U\ is the on-site attrac- 
tion leading to s-wave SC. We will measure all energies 
in units of t. At weak disorder this model has a super- 
conducting ground state for any density (n) ^ 1, whereas 
strong disorder (V ^> 1) localizes the fermions and de- 
stroys superconductivity. 




FIG. 1: Energy and temperature scales across the 
superconductor-insulator transition (SIT). The super- 
conducting T c (blue dots) decreases to zero at the critical 
disorder strength V c . The single-particle gap cjdos (black dia- 
monds), obtained from the DOS shown in Fig. 2, is large and 
finite in all states. The two-particle gap cj pa i r (red squares), 
obtained from the dynamical pair susceptibility shown in 
Fig. 3, is non-zero in the insulator but vanishes at the SIT. 
The dashed curves are guides to the eye; extracting critical 
exponents requires finite-size scaling beyond the scope of this 
paper. The statistical error bars in all the figures are dom- 
inated by disorder averaging and not from the QMC. These 
results were obtained at fixed attraction \U\ = 4 and aver- 
age density (n) ~ 0.87 on up to 100 disorder realizations on 
8x8 lattices, and uJdos were calculated at the lowest 

accessible T = 0.1. 

We use the determinantal QMC method [22], which 
is free of the fermion sign problem for the Hamiltonian 



(1), to compute various physical observables as functions 
of temperature and disorder. We work at \U\ = 4, so 
that the coherence length is within the system size, and 
work at a density (n) = 0.875. We have made exten- 
sive comparisons of the QMC results with self-consistent 
Bogoliubov-deGennes (BdG) calculations, which only 
take into account only the spatial amplitude variations; 
see supplementary material. These comparisons permit 
us to separate the effects of amplitude inhomogeneity and 
phase fluctuations. 

We compute frequency-dependent observables across 
the SIT for the first time. The single-particle DOS, 
LDOS and the pair susceptibility are obtained using the 
maximum entropy method (MEM) for analytic continua- 
tion [23, 24]. We have verified that these results obey var- 
ious sum rules to high precision, and that the MEM cor- 
rectly reproduces the low-energy structure of test spectra 
as shown in the supplementary material. What gives us 
confidence is that our central results on the single- and 
two-particle gaps can be equally well estimated directly 
from the exponential decay of the imaginary-time QMC 
data, without recourse to MEM. 

Phase diagram: In Fig. 1 we summarize our key re- 
sults for the disorder dependence of various temperature 
and energy scales. The Berezinskii-Kosterlitz-Thouless 
(BKT) critical temperature T c is estimated from the uni- 
versal jump in the superfluid density p s , calculated from 
the transverse current correlator [25]. We note that this 
procedure on finite systems provide an upper bound on 
the actual T c in the thermodynamic limit. As disorder 
strength V increases, T c falls and finally vanishes at the 
critical disorder V c , which defines the SIT. The single- 
particle energy gap uj^os remains non-zero across the SIT, 
whereas the two-particle energy gap cj pa i r is finite in the 
insulator and goes to zero at the transition. These gap 
scales are extracted from the DOS and the dynamical 
pair susceptibility, discussed below. Figure 1 can be in- 
terpreted as a phase diagram: T c is the superconducting 
transition temperature, co> pa i r is a crossover scale between 
the insulator and the quantum critical region, and c^dos 
is the pseudogap crossover scale described below. 

Single-particle spectra: We show in Fig. 2 the disor- 
der and temperature dependence of the DOS N(u). Pan- 
els (A,B) show the evolution with disorder at a very low 
temperature T = 0.1. The gap uJd os clearly remains fi- 
nite in both superconducting and insulating states. This 
counterintuitive prediction agrees qualitatively with BdG 
theory [13, 14]. Note, however, that the coherence peaks 
are predominantly present in the superconducting state 
and absent in the insulating state. Thus, it is the coher- 
ence peaks, and not the gap, that serve as spectral signa- 
tures of superconducting order. 

Figure 2(C,D) show the temperature evolution of N(uj) 
at weak disorder V < V c . Unlike BCS theory, the hard 
SC gap does not close with increasing T. Instead, the 
coherence peaks gradually disappear as the temperature 
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FIG. 2: The single-particle DOS N(uj) (upper panels) and representative spectra (lower panels), along three different cuts 
through the temperature-disorder plane. Panels (A,B): Disorder dependence of N(oj) at a fixed low temperature. A hard gap 
(black region) persists for all V above and below the SIT (V c ~ 1.6), but the coherence peaks (red) exist only in the SC state 
and not in the insulator. Panels (C,D): T-dependence of the N(cj) for the superconductor (V < V c ). The coherence peaks (red) 
visible in the SC state, vanish for T > T c w 0.14. A disorder-induced pseudogap (loss of low-energy spectral weight) persists 
well above T c . Panels (E,F): T-dependence of N(oj) for the insulator (V > V c ). The hard (insulating) gap at low T evolves 
into a pseudogap at higher T. No coherence peaks are observed at any T. 



increases across T c . Above T c , the gap gradually fills up, 
with a pseudogap persisting well above T c . 

The temperature evolution of N(uj) at strong disorder 
V > V c is shown in Fig. 2(E,F). Here the ground state 
is an insulator with a hard gap and little evidence for 
coherence peaks, and the pseudogap persists up to an 
even higher temperature. 

Two-particle spectra: What is the energy scale 
in the insulator that vanishes at the quantum criti- 
cal point? We propose that it is the gap for a two- 
particle excitation in the insulator. To access this 
gap, we examine the pair susceptibility P{oS) obtained 
by analytical continuation of the correlation function 
P{r) = £ R (T T F(R;T)Ft(R;0)) where F(R,r) = 

c iu( r ) c RT( r )' Thus P(t) is the amplitude for a pair cre- 
ated at a site R at r = to be found at the same site at 
a later time r. We find that in the insulating phase P(r) 
decays exponentially, which allows us to define uj pa ir, the 
characteristic energy scale for two-particle excitations. 

In Fig. 3 we show the imaginary part of the pair sus- 
ceptibility P"(uj) for three disorder strengths. At weak 
disorder P"(uj) is very large at low cj, whereas at strong 
disorder it has a clear two-peak structure with a char- 
acteristic energy scale cj pa i r . This dominant scale rep- 
resents the typical energy required to insert a pair into 
the system. We find that cj pa i r collapses to zero at the 



SIT because there is no cost for inserting a pair into a 
condensate. Note that there is the possibility of low en- 
ergy weight in P" (uj)/uj originating from rare regions but 
their spectral weight is small. 

Local probes: In Fig. 4 we track the behavior of var- 
ious quantities with increasing disorder strength]/. We 
show the LDOS iV(R, uj) at representative points, maps 
of the spatial variation of the density n(R), and the BdG 
pairing amplitude A(R) = (cr^cr^) (which cannot be 
computed in QMC). We see that the system becomes in- 
creasingly inhomogeneous with increasing disorder (mov- 
ing from left to right in Fig. 4). The SIT occurs due to 
loss of phase coherence between superconducting islands, 
seen as blue patches in the pairing amplitude A-map. 

We predict experimentally measurable signatures of 
the local density and pairing amplitude in the LDOS 
AT(R, uj). Let us focus on three representative sites Ri, 
R2, and R3. At moderate and strong disorder, Ri is 
located on a potential hill, with a low density n(Ri) w 
and a negligible pairing amplitude A(Ri) « 0. Thus the 
LDOS at Ri is highly asymmetric LDOS, with most of 
the spectral weight at uj > 0, for adding an electron. In 
contrast, R3 is in a potential well, with a high density 
ra (R3) « 2 and a negligible pairing amplitude A(R 3 ) « 0. 
Thus R3 also has a highly asymmetric LDOS, but most of 
the spectral weight is at uj < 0, for removing an electron. 
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FIG. 3: Imaginary part of the dynamical pair sus- 
ceptibility P"(u)/u at T = O.lt, averaged over 10 disorder 
realizations at three disorder strengths. Error bars represent 
variations between disorder realizations. For V <V C , there is 
a large peak at uj = 0, indicating zero energy cost to insert a 
pair into the SC. For V > V c , there is a gap-like structure with 
an energy scale cj pa ir, the typical energy required to insert a 
pair into the insulator which increases with V. 



Finally, R2 lies in a superconducting puddle close to 
half-filling, n(R,2) ~ 1, which permits particle mixing, 
and therefore a large pairing amplitude A(R 2 ). The 
LDOS at R2 is much more symmetrical, with large co- 
herence peaks that persists across the SIT and even in 
the insulating state. Note that all the LDOS curves have 
a clear gap. Thus, the more symmetrical coherence peaks 
in the LDOS, and not the local energy gap, are a clear 
experimental signature of local pairing. 

Conclusion: A detailed picture of the superconductor 
to insulator transition has emerged from our microscopic 
calculations. Even "homogeneous disorder" arising from 
an uncorrelated random potential leads to a highly inho- 
mogeneous state with superconducting puddles. These 
puddles shrink with increasing disorder, leading to en- 
hanced quantum phase fluctuations that destroy phase 
coherence at the SIT, resulting in an unusual insulator 



with a large single-particle gap but a much lower energy 
scale for pair excitations that vanishes at the SIT. In ad- 
dition, our results suggest that coherence peaks in the 
DOS are destroyed by thermal fluctuations for T > T c , 
but a pseudogap persists well above T c . This disorder- 
induced pseudogap is a robust consequence of the small 
superfluid stiffness [26] and would exist near the SIT even 
in a highly disordered weak-coupling system. 

Our goal was to model the disorder-driven SIT in s- 
wave SC films, but aspects of our results bear a striking 
resemblance to the completely different - and not well 
understood - problem of the pseudogap in underdoped 
cuprates. Features like the loss of low-energy spectral 
weight persisting across thermal or quantum phase tran- 
sitions, even as coherence peaks are destroyed, may well 
be common to all systems where the small superfluid stiff- 
ness drives the loss of phase coherence. The cuprate pseu- 
dogap is driven by the proximity to the Mott insulator 
and competing order parameters, with disorder proba- 
bly playing a secondary role, unlike the disorder-induced 
pseudogap near the SIT discussed in this paper. 
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FIG. 4: Local density of states (LDOS) iV(R, lj), density n(R), and BdG pairing amplitude A(R) as a function of disorder 
strength for a montage of nine disorder realizations of 8 x 8 lattices. Panels A, B, C correspond to V = 0.1, 1, 2 respectively. 
The LDOS is plotted at three representative sites Ri. At moderate and strong disorder, site Ri is on a high potential hill that 
is nearly empty, and R3 is in a deep valley that is almost doubly occupied. This leads to the characteristic asymmetries in the 
LDOS in the center and right columns for Ri and R3. The small local pairing amplitude A(R) at these two sites is reflected 
in the absence of coherence peaks in the LDOS. In contrast, site R2 has a density closer to half- filling, leading to a significant 
local pairing amplitude, a much more symmetrical LDOS, and coherence peaks that persist even at strong disorder. 
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SUPPLEMENTAL ONLINE MATERIAL 



PACS numbers: 



In this supplement we provide details of the determi- 
nantal QMC simulations, comparison between QMC and 
inhomogeneous Bogoliubov de-Gennes (BdG) mean field 
theory, and the analytic continuation procedure for ex- 
tracting real frequency information from imaginary time 
QMC data. 

Determinantal QMC: We use the determinan- 
tal Quantum Monte Carlo (QMC) algorithm [1, 2] to 
calculate the quantities discussed in the paper, in- 
cluding the imaginary-time Green function G(R; r) = 

^^ c Ra( r ) c Rcr(0)) an d pairing correlation function 
P(t) = En(TrF(K;r)F^K;0)) where F(R,r) = 
c R|( r ) c Rt( r )* We present results for 8 x 8 square lat- 
tices with periodic boundary conditions. The lattice size 
is dictated by the need for very accurate QMC data re- 
quired for analytic continuation. 

For a given set of parameters, the simulations are equi- 
librated for up to 4 x 10 5 Monte Carlo steps. The final 
averages for a single disorder realization are taken over 
2 x 10 5 steps for static quantities and over 4 x 10 6 for 
dynamical quantities. We further average over 10 disor- 
der realizations for a given disorder strength. The re- 
sulting maximum absolute errors are SG(r) ~ 10 -4 and 
SP(r) ~ 10- 2 . 

Comparisons of QMC with BdG: In Fig. 1 we 
show a comparison of the local density n(R) obtained 
using QMC and self-consistent BdG, including inhomo- 
geneous Hartree shifts, for one disorder pattern at dif- 
ferent disorder strengths. The close agreement indicates 
that phase fluctuations, not included in BdG, have very 
little effect on n(R). On the other hand, the superfluid 
stiffness and spectral properties at finite temperatures 
and large disorder are greatly affected by thermal and 
quantum phase fluctuations. 

The local density is directly related to the occupied 
and unoccupied part of the LDOS (see Fig. 4 of the pa- 
per) via the sum rules: 2 dojf(oj)N(H^uj) = n(R) 
and 2 du[l - f(u)]N(R, u) = 1 - n(R), where f(u) 
is the Fermi function and the factor of 2 comes from spin 
degeneracy. We have tested these sum rules for the calcu- 
lated LDOS and find excellent agreement. Further sum 
rule tests are described below. 

Analytic continuation: We use the maximum en- 
tropy method (MEM) to extract the local density of 
states iV(R, uo) and the pair spectrum P"(uj) from the 
imaginary-time Green function G(R; r) and pairing cor- 
relation function P(r) respectively. The MEM for an- 
alytic continuation [3] essentially inverts the Laplace 
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FIG. 1: Comparison of the local density n(R) obtained using 
QMC (top panels) and self-consistent BdG (middle panels) 
for one disorder pattern at three different disorder strengths. 
The densities are very similar as seen from the histograms of 
their differences in the lowest panels, indicating that phase 
fluctuations have very little effect on n(R). 
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The average DOS N{uS) is obtained from analytic contin- 
uation of G(R, r). 

We have performed extensive tests using known model 
spectra as follows: (i) choose a test spectrum N(uj); (ii) 
perform a Laplace transform to obtain the imaginary- 
time Green function G(r); (hi) add random numbers 
5G(r) drawn independently from a normal distribution 
of width SG = 10 -4 , in order to simulate Monte Carlo 
statistical error; and finally (iv) feed the resulting noisy 
data, Gdata(^), into our MEM routine. This procedure is 
illustrated in Fig. 2. We have concluded that the MEM 
is adequate for extracting the low-energy features of the 
spectrum, particularly the gap. 

We emphasize that two of our crucial observations re- 
garding the single particle gap scale in N(uj) and the 
two-particle gap scale in P"(uj)/uj can be directly esti- 
mated, from the exponential decay at large r of G(r) 
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FIG. 2: Demonstration of analytic continuation using the maximum entropy method (MEM). A test spectrum (blue) is chosen, 
noise is added in the r-domain, and the MEM is used to reconstruct the spectrum (red). The reconstructed Green function 
Gmem(t) fits the "data" Gdata(r) to within its error bars. The reconstructed spectrum Nmem(w) agrees well with the test 
spectrum N(oj) at low energies, reproducing the correct gap structure, although there are deviations at higher energies. Indeed, 
the magnitude of the gap can be estimated by examining the exponential decay constant of G(r), even without using the MEM. 



and P(t) calculated from the QMC data without any 
MEM analytic continuation. 

We have also made extensive sum-rule checks for the 
spectra obtained from MEM. We define the rath 

frequency moment of the local density of states at posi- 
tion R, as 
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These moments satisfy the following sum rules, which can 
be derived rigorously by extending the analysis in Ref. [4] 
to a disordered system: 
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where z = 4 is the coordination number. Differentiating 
Eq. (1) shows that the moments are also related to the 
values and derivatives of the Green function at r = and 
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QMC simulations produce very accurate results for 
Gr(0) and Gr(/3), with absolute errors of about 10 -5 . 
The MEM analytic continuation procedure fits these data 
points to within the error bars. We have verified that the 
MEM LDOS satisfies the moment sum rules with a frac- 
tional error of less than 0.001% for M®\ less than 1% 
for and about 1% for . In contrast, the LDOS 

obtained from BdG calculations are found to violate the 
sum rules. 
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